Free energy surface of ST2 water near the liquid-liquid phase transition 
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We carry out umbrella sampling Monte Carlo simulations to evaluate the free energy surface of 
the ST2 model of water as a function two order parameters, the density and a bond-orientational 
order parameter. We approximate the long-range electrostatic interactions of the ST2 model using 
the reaction-field method. We focus on state points in the vicinity of the liquid-liquid critical point 
proposed for this model in earlier work. At temperatures below the predicted critical temperature 
we find two basins in the free energy surface, both of which have liquid- like bond orientational order, 
but differing in density. The pressure and temperature dependence of the shape of the free energy 
surface is consistent with the assignment of these two basins to the distinct low density and high 
density liquid phases previously predicted to occur in ST2 water. 



I. INTRODUCTION 

In 1992, the results of a computer simulation study of 
the ST2 model [l[ of water were used to propose that 
a liquid-liquid phase transition (LLPT) occurs in super- 
cooled water [2J. Below the critical temperature T c for 
the proposed LLPT, two distinct phases of water, the 
low density liquid (LDL) and high density liquid (HDL) 
phases are separated by a first-order phase transition. 
The predicted phase diagram for the ST2 model in the 
plane of temperature T and pressure P in the vicinity of 
the critical point is shown in Fig. [TJ 

An appealing feature of the LLPT proposal is that it 
simultaneously accounts for (a) the unusual thermody- 
namic behavior of liquid water in the supercooled region, 
and (b) the occurrence of two distinct forms of amor- 
phous solid water in the glassy regime [!, H| . Evidence 
for a LLPT has been reported in a number of simulation 
studies of water and water- like systems; see e.g. [oT-frij]. 
Experimentally, a LLPT has yet to be decisively con- 
firmed in supercooled water, and efforts to resolve this 
question in the laboratory continue (l2Tll5j ]. The pre- 
dicted location of the critical point in the supercooled 
regime is challenging to study in experiments because of 
rapid ice crystallization. In simulations, this problem is 
avoided when the liquid can be studied on a time scale 
that is long relative to the liquid-state relaxation time, 
but short compared to crystal nucleation times. 

Recently, Limmer and Chandler [l6| have challenged 
the LLPT hypothesis. Using umbrella sampling Monte 
Carlo (MC) simulations of two water models (mW [17| 
and ST2 water), Ref. [l6[ presents results for the free en- 
ergy surface F(p,Qe) of the liquid as a function of two 
order parameters, the density p, and a bond-orientational 
order parameter Qe- Qe is a bulk order parameter used 
to distinguish crystalline configurations from liquid or 
amorphous solid states of a system. Values of Qe ap- 
proaching zero correspond to disordered states, while 
larger values of Qe indicate greater degrees of crystalline 



order. The detailed definition of Qe is given in Eqs. 1-3 
of Ref. [l6j], and is based on an analysis of the orienta- 
tion of local molecular environments (i.e. a molecule and 
its nearest neighbors) in terms of spherical harmonics, 
as originally proposed by Steinhardt, et al. (is) . In the 
present work, we use the same definition of Qe as given 
in Ref. 0. 

It has long been appreciated that the density of the 
proposed LDL and HDL phases must be different. The 
innovation of Ref. [l6| is that by examining the depen- 
dence of F(p,Qe) on Qe, Limmer and Chandler address 
the relationship of the metastable liquid phase to the or- 
dered crystalline ice phases. If a LLPT transition occurs 
in a simulation model, then under appropriate conditions 
of T and P, two distinct free energy basins should be ob- 
served in F{p,Qe) in the \cw-Qe (i.e. liquid-like) regime. 
For both the mW and ST2 water models, Ref. [16| re- 
ports that only one liquid-like free energy basin is found 
in F(p, Qe), including, in the case of ST2 water, at con- 
ditions below the proposed critical temperature of the 
LLPT. Limmer and Chandler conclude that phenomena 
previously interpreted as evidence for a LLPT are in fact 
due to the liquid-to-crystal phase transition. 

Since the publication of Ref. Liu et al. have re- 
ported on their own evaluation of the free energy surface 
F(p,Qe) found from umbrella sampling MC simulations 
of ST2 water [l9j]. Althou gh t hey employ methods sim- 
ilar to those used in Ref. jl6l|. Liu et al. report a very 
different result: the observation of two distinct liquid 
free energy basins in F(p, Qe), with properties consistent 
with the LLPT hypothesis. The results of Ref. [lj| are 
also consistent with an earlier study by the same group 
reporting the free energy of ST2 water as a function of p 
only [H- 

The precise reasons for the difference between the re- 
sults of Refs. and [H for F(p,Q 6 ) remain unclear. 
Among the differences in the approaches used in these 
two works, we note two. First, Limmer and Chandler 
present results for F(p, Qe) at various pressures as deter- 
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FIG. 1: Phase behavior of the ST2-RF model predicted from 
previous work using N = 1728 molecular dynamics simu- 
lations. Shown are the estimated locations of the critical 
point (green circle) and the LDL-HDL coexistence line (green 
line) [24|. Note that the error bars associated with the crit- 
ical point also apply to the coexistence line. Estimates for 
the HDL spinodal (down-triangles) and LDL spinodal (up- 
triangles) are also shown [f|. Red circles locate the state 
points at which we carry out series K, L, and M of the present 
work. 



mined at one temperature, T — 235 K, which is below 
but within error of the estimated critical temperature 
T c = 237 ± 4 K for the ST2 model when studied with 
Ewald summations [20]. Working this close to T c may 
make it difficult to discern distinct liquid basins in the 
free energy surface within the statistical error. Liu et 
al. report results for a range of temperatures below T c , 
from T = 224 to 235 K, and show that the distinction be- 
tween the two liquid basins that they observe in F(p, Q§) 
becomes greater as T decreases below T c . 

Second, in both Refs. [l|| and [Hj], the method of 
Ewald summation is used to approximate the long-range 
contributions to the electrostatic potential energy of the 
ST2 system. However, Liu et al. report that their 
Ewald summation method employs vacuum boundary 
conditions, whereas Limmer and Chandler use conduct- 
ing boundary conditions. Liu et al. note some significant 
sensitivity in the behavior of their system as a function 
of these boundary condition choices. If and how these 
boundary conditions might affect the qualitative shape 
of the F(p,Q 6 ) surface is incompletely understood. 

In light of the conflicting results of Refs. and (l9j . 
we present here a new evaluation of the free energy sur- 
face F(p,Q 6 ) of ST2 water. In order to expand our 
understanding of the role of long-range interactions, we 
use a different approach to account for the electrostatic 
energy namely the reaction field method [2lj . Indeed, 
many of the previous studies of ST2 water that relate to 
the LLPT hypothesis were conducted using the reaction 
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FIG. 2: Time dependence of the collective intermediate scat- 
tering function f(t) for runs with various values of p* in series 
(a) K and (b) L. Each curve is an average over 10 runs. 



field method [H, 0, [23| , including the work in which 
the occurrence of a LLPT was first proposed [2J . Further- 
more, a recent umbrella sampling MC study of the ST2 
model, using the reaction field method, showed that the 
shape of the free energy as a function of p was consistent 
with the LLPT hypothesis [22j. An explicit examina- 
tion of the F(p, Qq) surface for the ST2 model with a 
reaction field treatment of the electrostatics is therefore 
warranted. In addition, we also study a range of tem- 
peratures and pressures in the vicinity of the proposed 
critical point, to examine their influence on the F(p. Qq) 
surface. 
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FIG. 3: Collective intermediate scattering function f(t) for the lowest value of p* from each series: (a) K, p* — 0.93 g/cm 3 ; 
(b) L, p* = 0.95 g/cm 3 ; and (c) M, p* = 0.95 g/cm 3 . These are the most slowly relaxing runs used in our analysis. The black 
lines give f(t) for each of the 10 independent runs conducted at the same values of (T, P, p* ,Q%), and the thick red line is their 
average. 



II. ST2 MODEL 

We study the ST2 model of water proposed by Still- 
inger and Rahman [![. The ST2 pair potential is a 
sum of a Lennard- Jones (LJ) interaction (centered on 
the O atom) , and electrostatic interactions involving four 
tetrahedrally positioned charges. Our model parameters 
for the geometry and pair interactions of the ST2 water 
molecule are the same as those given in Ref. The 
potential energy U of our system is given by, 



C/ = C/ e + (7lj + A ( 7 L j, 



(1) 



where Us, and E/lj are the respective electrostatic and LJ 
contributions. In our simulations, the LJ interaction is 
sharply cut off when the 0-0 distance r exceeds R c = 
0.78 nm, and the contribution from longer ranged LJ 
interactions is approximated by, 



AC/lj = 



8ne(T 6 p„N 
3M ' 



(2) 



as described in the Appendix of Ref. In Eq. [51 N 
is the number of molecules, p n is the number density of 
molecules, and e and a are the respective energy and size 
parameters of the LJ potential. 

To evaluate /7e, the electrostatic contributions to the 
potential energy, we adopt the treatment used in the 
study of ST2 water by Steinhauser; see Eqs. 5 and 6 of 
Ref. [2lJ]. In this approach, the electrostatic interactions 
of the ST2 model are evaluated directly up to r — R c us- 
ing the original form given in Ref. including the use 
of a "switching function" to preclude a divergence of the 



energy due to charge overlaps. The contribution of elec- 
trostatic interactions beyond R c is then approximated 
using the reaction field method, in which the liquid be- 
yond R c is treated as a polarizable dielectric continuum. 
As in Ref. [2l| , we assume that the dielectric constant of 
the continuum liquid is £r , = 00. To avoid a sharp discon- 
tinuity in the electrostatic interactions at i? c , a tapering 
function (described in Ref. [2l|) is used to smoothly re- 
duce the electrostatic interaction between two molecules 
(both direct and reaction field contributions) to zero over 
the interval 0.95i? c < r < R c . 

The evaluation of the pair interactions as described 
above is the same procedure that was used in a number 
of previous studies 0, d, IE HH, HH H|| . For the remainder 
of this paper, we will refer to the reaction field version 
of ST2 described above as ST2-RF, to emphasize the dif- 
ference between the present study and those works that 
have studied the ST2 model using an Ewald treatment 
of the electrostatics flU fl9ll20ll . 



III. SIMULATION METHODS 

Our aim is to evaluate the free energy surface F(p, Qe) 
for the ST2-RF model in the vicinity of the predicted 
LLPT for this model. To define F(p,Q e ), let p{p,Q e ) 
be proportional to the equilibrium probability for a mi- 
crostate of the system at fixed values of N, T, and P to 
have order parameter values p and Qq. The conditional 
Gibbs free energy F(p, Q 6 ) is then defined by, 



F(p,Q 6 



-kT\np(p,Q 6 ) + F Q 



(3) 
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FIG. 4: Comparion of r a (filled symbols) and r run (open sym- 
bols) as a function of p* for series K (circles), L (squares), and 
M (diamonds). Values of r a < 10 4 MCS are not shown. 



where Fq is an (irrelevant) constant related to the nor- 
malization of p, and k is Boltzmann's constant [l6j . 

We also define the "contraction" of F with respect to 
Qe as, 



HP) 



-fcTln 



<?;: 



dQ 6 exp[-/3F(p, Q 6 )} 



(4) 



where /3 = 1/kT [T(|. F(p) represents the free energy as 
a function of p that would be found from an ensemble of 
states in which Q 6 is free to vary between zero and Qjf ax . 
In this work, we are concerned with the liquid-like range 
of Qq. As shown below, we find that setting Q™ ax = 0.09 
is sufficient to characterize F{p) for the liquid-like basins 
of the free energy surface. 

Following the approach of Ref. [llj], we use umbrella 
sampling MC simulations to evaluate F(p, Qe)- We carry 
out MC simulations in the constant- (N, P, T) ensemble, 
and to implement umbrella sampling, we add a biasing 
potential, 



U B =ki(p- p* 



fe(Q 6 -Qe) 2 



(5) 



to the system potential energy U in Eq. [T] The effect of 
Ub is to constrain a given simulation to sample configu- 
rations in the vicinity of chosen values of the order pa- 
rameters p = p* and Qe = Qq. In all our simulations, we 
fix TV = 216, fci = 1000/cT (cm 3 /g) 2 , and k 2 = 2000/cT. 

Trial configurations for each Monte Carlo step (MCS) 
are generated as follows: First, we carry out a mini- 
trajectory of 10 unbiased (i.e. Ub = 0) constant- 
(N,P,T) MC moves, in which each move consists (on 
average) of N — 1 attempted rototranslational moves, and 
one attempted change of the system volume. The maxi- 
mum size of the attempted rototranslational and volume 
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FIG. 5: Contracted free energy F(p) at T = 240 K, at two 
different pressures. The filled circles are obtained by analyz- 
ing series K and M at P = 215 MPa, the pressure at which 
these series are conducted; the open circles are obtained by 
reweighting these results to P — 200 MPa. The filled squares 
are obtained by analyzing series L at P = 200 MPa, the pres- 
sure at which this series is conducted; the open squares are 
obtained by reweighting these results to P = 215 MPa. 



changes are chosen to give MC acceptance ratios in the 
range 25-40%, depending on the thermodynamic condi- 
tions. Next, the change in the biasing potential Ub is 
evaluated for the trial configuration resulting from the 
mini-trajectory, relative to the system configuration at 
the beginning of the mini-trajectory, to determine the 
acceptance or rejection of the trial configuration. This 
completes one MCS, and the procedure is then repeated. 

In order to identify the T-P state points at which 
to conduct our runs, we use the location of the LLPT 
reported in previous work. Fig. [1] shows the estimates 
for the critical point and coexistence line obtained from 
N = 1728 molecular dynamics simulations of the ST2-RF 
model. Of particular importance are the locations of the 
spinodal lines for the LDL and HDL phases. These spin- 
odal lines demarcate the stability limits for each phase. 
Consequently, if liquid-liquid coexistence does indeed oc- 
cur in the ST2-RF model, the F(p,Qe) surface will si- 
multaneously exhibit two distinct liquid basins only for 
state points lying in the region between the HDL and 
LDL spinodals. It is in this region of states that we fo- 
cus our simulations. To carry out our runs, we select 
pressures that lie between or near to the HDL and LDL 
spinodals, and a temperature (T = 240 K) that is 7 K 
below the estimated critical temperature of T c — 247 ± 3 
for the ST2-RF model [24 1. 

We carry out three distinct series of runs. In the fol- 
lowing, "series K" denotes the set of runs conducted at 
T = 240 K,P= 215 MPa, Q* e = 0.05, and equally 
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FIG. 6: Contracted free energy F(p) at T = 240 K and P = 
204.5 MPa, obtained by combining all results from series K, 
L, and M. 



spaced values of p* from 0.93 to 1.15 g/cm 3 , separated 
by 0.01 g/cm 3 . "Series L" denotes runs conducted at 
T = 240 K, P = 200 MPa, Q* 6 = 0.05, and equally 
spaced values of p* from 0.95 to 1.15 g/cm 3 , separated by 
0.01 g/cm 3 . "Series M" denotes runs conducted at T = 
240 K,P = 215 MPa, Q* 6 = 0.09, and p* = 0.95 g/cm 3 . 
The state points in the T-P plane corresponding to se- 
ries K, L, and M are identified in Fig. [T] For all distinct 
choices of (T,P,p*,Qg) in the above series, we conduct 
10 separate runs, each initiated from independent start- 
ing configurations. The results presented here are thus 
based on an analysis of 450 independent runs. 

All our runs are carried out for between 5 x 10 6 to 
5 x 10 7 MCS. Using the second half of each run, we com- 
pute /(t), the collective intermediate scattering function 
as a function of time t. We evaluate fit) at the lowest- 
wavenumber peak in the static structure factor for the 
O atoms, i.e. the so-called first sharp diffraction peak of 
molecular tetrahedral networks. As shown in Figs. [2] and 
O in all cases f(t) decays to zero on a time scale which 
is short compared to the lengths of our runs. Hence the 
system behaviour is consistent with liquid-like relaxation 
under all conditions simulated in this study. After aver- 
aging f(t) over the 10 runs at each choice of (T, P, p* , Q$), 
we estimate the alpha-relaxation time Tq clS the time at 
which f(t) = er l . As shown in Fig. 21 in all cases we 
find r a < 2 x 10 5 MCS. To account for equilibration, we 
then discard the results for t < r e of each run, where 
r e = 20t q or 10 4 MCS, whichever is larger. The result- 
ing length r run of each production run that is used in 
our analysis is shown in Fig. |4l compared to the corre- 
sponding value of r a . In terms of r a , the lengths of our 
production runs range between 175r Q and 4400r Q . As 
shown in Fig. [3J even our most slowly relaxing individual 



simulations are run for a time that is at least two orders 
of magnitude longer than the corresponding value of r a . 

To estimate F(p,Q 6 ), F(p), and the associated error, 
we use the multistage Bennet acceptance ratio (MBAR) 
method [26| . The MBAR method takes as input the time 
series of the order parameters (/? and Qe) and the sys- 
tem potential energy U, reweights the statistics obtained 
from each run to remove the effect of the biasing po- 
tential, and produces an optimal estimate of the desired 
free energy function at a specified value of T and P. The 
MBAR method also facilitates reweighting the configura- 
tions sampled during our runs with respect to T and/or 
P, allowing the statistics from different state points to be 
combined to produce an estimate of F(p,Q 6 ) or F(p) at 
T-P state points that lie near to the conditions at which 
we carry out our simulations. 

For the purpose of estimating the free energy and its 
error using MBAR, we wish to consider only those config- 
urations from our runs that are statistically independent. 
We assume that statistically independent configurations 
are separated by r a or 10 4 MCS, whichever is larger. All 
other configurations are ignored in our analysis. Note 
that in all our plots the indicated error is the error with 
respect to the minimum value of the estimated free en- 
ergy, which in most cases is arbitrarily set to zero. Also, 
all error bars reported here represent one standard devi- 
ation of error. 



IV. RESULTS 

First, we compare the results obtained for F(p) at 
the two state points directly simulated in our runs. Se- 
ries K and M are both conducted at T — 240 K and 
P = 215 MPa, while series L is conducted at T = 240 K 
and P — 200 MPa. The results for F(p) obtained using 
only series K and M, and that obtained using only series 
L are compared in Fig. [5] The shapes of both curves sug- 
gest the existence of two distinct free energy minima sep- 
arated by an interval of thermodynamic instability with 
respect to p, as indicated by concave-down curvature of 
F(p). One minimum is centred near 0.9 g/cm 3 and the 
other near 1.05 g/cm 3 . 

To check that the statistics we have gathered in se- 
ries K and M are consistent with the results obtained 
from series L (and vice versa), we also show in Fig. O 
the result for F(p) found by reweighting our data from 
series K and M to P — 200 MPa, and the result found 
by reweighting our data from series L to P = 215 MPa. 
The reweighted results are in good agreement with the 
unreweighted curves, confirming that both data sets have 
independently converged to equilibrium. In the remain- 
der of this paper, all results shown for F(p) and F(p, Qq) 
are therefore obtained by combining the statistics from 
all three simulations series, K, L, and M. 

In Fig. [U we show F(p) at T = 240 K and P — 
204.5 MPa, a pressure intermediate between those shown 
in Fig. [51 At this state point, F(p) clearly displays two 
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FIG. 7: Contour plots of F(p,Q e ) at T = 240 K for (a) P = 195 MPa, (b) P = 204.5 MPa, and (c) P = 230 MPa. To 
evaluate these surfaces, we have coarse-grained the plane of p and Qe into rectangular cells of dimensions Ap = 0.02 g/cm 3 and 
AQa = 0.01. Data from all series (K, L, and M) are combined and analyzed to obtain these plots. For each panel, contours 
are separated by 0.5fcT, and the error in F is 0.5/cT or less. The lowest lying values of F in each plot are labelled LDL and/or 
HDL. 



distinct free energy minima, separated by a free energy 
barrier of approximately lkT, a typical value when T is 
close to T c . 

We next analyze the behavior of the free energy sur- 
face F(p, Qe). Fig. shows contour plots of F(p, Qq) at 
T = 240 K for three pressures from P = 195 to 230 MPa. 
The F(p, Qq) surface at P = 204.5 MPa simultaneously 
displays two free energy basins, each corresponding to 
a distinct metastable thermodynamic phase. The min- 
ima of both basins are located at liquid-like values of 
Qq in the range 0.05-0.065. The shape of both basins 
shows that the phases they represent are locally stable 
with respect to fluctuations in both p and Qq. The sta- 
bility of both phases with respect to Qe, highlighted in 
Fig. [HI shows that neither free energy basin is connected 
via a monatonic "downhill" path to any of the free en- 
ergy basins associated with the various phases of crys- 
talline ice, which are expected to occur at much higher- 
values of Qq ~ 0.5. The properties of the phases associ- 
ated with the two basins shown in Fig.[7[b) are therefore 
consistent with two distinct liquids, the LDL and HDL 
phases, predicted to occur in the ST2-RF model in earlier 
work @, 0, HI ■ 

If the two basins shown in Fig. Efb) are consistent with 
a LLPT between LDL and HDL phases, then increasing 
the pressure at constant T should cause the LDL basin 
to disappear, and decreasing the pressure should cause 
the HDL basin to disappear, as both phases reach the 
respective spinodal limits that bracket the coexistence 
curve (see Fig.Q}. This is illustrated in Fig.JTJa) and (c). 



At P = 195 MPa, only the LDL basin remains, while at 
P = 230 MPa only the HDL basin is observed. 

We note that at T = 240 K, the pressure range found 
here that corresponds to the region between the HDL and 
LDL spinodals appears to be shifted downward by about 
10 MPa relative to the thermodynamic features shown 
in Fig. [T] However, this difference is less than the error 
associated with the results in Fig. Q] for the location of 
the critical point and coexistence line. Considering that 
the features in Fig. Q] are based on an extrapolation of 
equation-of-state data from N = 1728 molecular dynam- 
ics simulations @, HiJ , and considering the possibility of 
differences due to finite-size effects when comparing with 
our N — 216 results, the agreement between the behav- 
ior observed here and that predicted in Fig. [T] is quite 
satisfactory. 

Finally, in Fig. |H]we show the evolution of F{p) along 
a path in the T-P plane that approaches the vicinity of 
the predicted critical point in ST2-RF. Consistent with 
the occurrence of a line of first-order phase transitions 
terminating in a critical point, the two basins in F(p) 
are separated by a higher free energy barrier at lower T, 
which decreases in height, and then disappears, on ap- 
proach to the critical point. Fig. [§] also confirms that 
the density of the HDL phase varies significantly with T, 
whereas that of the LDL phase is comparatively insensi- 
tive to changes in T . This observation is consistent with 
previous studies of the free energy of ST2 water that have 
observed distinct HDL and LDL basins [HIM El- 
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FIG. 8: Slices through the free energy surface F(p,Qa) for 
T — 240 K and P = 204.5 MPa [shown in Fig. □»] as 
a function of Qr,, passing through the minima of the LDL 
basin at p — 0.90 g/cm 3 (circles), and the HDL basin at 
p — 1.04 g/cm 3 (squares). 



V. DISCUSSION 

In summary, for the ST2-RF model, we find two dis- 
tinct basins in the free energy surface F(p, Qe), differing 
in density, but both occurring at low values of Qq, assur- 
ing that they correspond to disordered thermodynamic 
phases. Furthermore, our results for the structural re- 
laxation times demonstrate that both basins correspond 
to equilibrated metastable liquid phases. These observa- 
tions, and the dependence of the shape and position of 
the basins as a function of T and P are entirely consistent 
with the occurrence of a LLPT in the ST2-RF model of 
water, as described in previous work 0, H, H, [H, HH . Our 
results are also consistent with those of Liu and coworkers 
for the ST2 model using an Ewald treatment of the elec- 
trostatics [ljl H(| . Our results are qualitatively different 
from the behavior of the ST2 system reported by Limmcr 
and Chandler [161 ] , and also are not consistent with their 
proposal that all the behavior previously ascribed to a 
LLPT in water-like models is in fact associated with the 
liquid-to-crystal transition. 

We note that Limmer and Chandler have argued that 
the observation of two liquid basins in F(p, Qq) could 
arise as an artifact of restrictin g th e sampling to low 
values of Qe; see Fig. 9 of Ref. [27| and the accompa- 
nying discussion. Limmer and Chandler note that their 
own data for F(p) "exhibits an inflection or slight mini- 
mum" for low values of Q™ ax , but that this shoulder in 
the curve merges into the minimum associated with the 
crystal basin for larger values of Qjf ax . From this behav- 
ior they conclude that although the shoulder observed 
for small Q™ ax "could be confused with a second liquid 
basin," it is in fact "due to the barrier separating liquid 
from crystal." (27j . 




FIG. 9: Contracted free energy F(p) at several state points 
approaching the liquid-liquid critical point. From bottom to 
top, the state points are: T = 230 K and P = 245 MPa; 
T — 235 K and P = 225 MPa; T = 240 K and P = 
204.5 MPa; T = 245 K and P = 184 MPa; and T = 250 K and 
P — 164 MPa. Each curve has been shifted by an arbitrary 
constant to facilitate comparison. The combined data from 
all series (K, L, and M) are analyzed to obtain each curve. 



We disagree with this interpretation of the data. We 
refer the reader to the bottom right-hand panel of Fig. 9 
of Ref. [27|, which shows the free energy surface upon 
which the above analysis of Limmer and Chandler is 
based. In this free energy surface, the liquid basin is 
clearly distinct from the crystal basin, in the sense that 
any path connecting the minima of these two basins must 
pass over a barrier of at least 23fcT. The shoulder in the 
free energy surface noted by Limmer and Chandler oc- 
curs deep inside the liquid basin (near p — 0.92 g/cm 3 
and Qq — 0.08), and is well separated from the barrier 
that defines the boundary between the liquid and crystal 
basins (near Qq = 0.27). Hence, any path leading from 
the shoulder to the crystal basin must go "uphill" in free 
energy at some point along the path. The shoulder thus 
cannot be understood as an extension of the crystal basin 
into the low-Q6 regime. When F(p) is plotted for differ- 
ent Q™ ax , the shoulder and the crystal minimum become 
superimposed on one another because they happen to oc- 
cur at similar densities; however, this is not a basis for 
concluding that these two features must be associated 
with the same (crystalline) free energy basin. 

To conclude, we emphasize that not all models of wa- 
ter exhibit a LLPT. For example, the mW model seems 
to be a case in which a LLPT, which otherwise might 
occur, becomes unobservable due to the loss of stability 
of the supercooled liquid with respect to crystal nucle- 
ation [IE |H 111. Whether or not a LLPT occurs in a 
given water model, and indeed in real water itself, may 
depend sensitively on the details of the intermolecular 
interaction. For real water, it remains for experiments to 
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determine conclusively if a LLPT can be observed, for ex- 
ample, by manipulating the rate of ice crystallization in 
the supercooled regime by exploiting nano- confinement 
or external fields. Nonetheless, the results presented here 
provide clear evidence that a LLPT does occur in simu- 
lations of the ST2-RF model of water, and confirm the 
conclusions drawn in previous studies of this model re- 
garding the existence of a LLPT. 
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